Orbital magnetization in periodic insulators 
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Working in the Wannier representation, we derive an expression for the orbital magnetization of 
a periodic insulator. The magnetization is shown to be comprised of two contributions, an obvious 
one associated with the internal circulation of bulk-like Wannier functions in the interior, and an 
unexpected one arising from net currents carried by Wannier functions near the surface. Each 
contribution can be expressed as a bulk property in terms of Bloch functions in a gauge-invariant 
way. Our expression is verified by comparing numerical tight-binding calculations for finite and 
periodic samples. 
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PACS numbers: 75.10.-b, 75.10.Lp, 73.20.At, 73.43.-f 

Recent years have seen a surge of interest in issues 
of charge and spin transport in magnetic materials and 
nanostructures, notably the development of a theory of 
the intrinsic anomalous Hall conductivity and some con- 
troversies surrounding the spin-Hall effect [1]. In this 
context it is quite surprising that the theory of orbital 
magnetization, essential for any proper description of 
magnetism, has remained in a primitive state. Linear- 
response methods allow calculations of magnetization 
changes [2-5] , but not of the magnetization itself. 

Hirst [6] has emphasized that a knowledge of the bulk 
local current density J(r) is insufficient, even in princi- 
ple, to determine the macroscopic orbital magnetization 
M, just as the density p(r) cannot be used to determine 
the electric polarization P. Thus, the theory of M to- 
day is in a condition very similar to that of P in the 
early 1990s, when the problem of computing finite po- 
larization changes was solved by the introduction of the 
Berry-phase theory [7, 8]. The essential difficulty, that 
the matrix elements of the position operator r are not 
well-defined in the Bloch representation, could be over- 
come by reformulating the problem in the Wannier rep- 
resentation. Because Wannier functions (WFs) are expo- 
nentially localized in an insulator, matrix elements of r 
between WFs are indeed well-defined. 

Here we show that it is possible to formulate a corre- 
sponding theory of the orbital magnetization for an insu- 
lator with broken time-reversal symmetry. The problem 
is analogous, with the circulation operator r x v now 
being ill-defined in the Bloch representation. Working 
instead in the Wannier representation, we write the or- 
bital magnetization as a gauge-invariant Brillouin-zone 
integral over occupied Bloch functions. It contains two 
terms, the first of which describes the internal circulation 
of bulk-like WFs [9]. The second is much more subtle, 
arising only from surface WFs and reflecting the fact that 
the information about surface currents needed to define 
the macroscopic magnetization is actually contained in 
the bulk bandstructure. The resulting formula is consis- 



tent with a recent semiclassical argument [10] and can 
easily be implemented in existing first-principles codes. 

For our derivation, wc restrict ourselves to the case 
of an insulator described by a one-particle Hamiltonian 
with broken time-reversal symmetry. While the restric- 
tion to insulators is essential for the theory of polariza- 
tion, we suspect that it is less so here, so that future 
generalizations to metals are not ruled out. We also re- 
quire a vanishing macroscopic magnetic field (or, more 
generally, an integer number of flux quanta per unit cell) 
so that the Bloch wavevector k remains a good quantum 
number. We have in mind cases in which a staggered 
magnetic field averages to zero over the unit cell, or in 
which the time-reversal breaking comes about through 
spin-orbit coupling to a background of ordered local mo- 
ments [11-15]. For simplicity we work with spinless elec- 
trons (the generalization to the spin-unrestricted case be- 
ing straightforward) and furthermore restrict ourselves to 
zero-Chern- number insulators [11, 12]. 

Let us consider a finite sample representing a fragment 
of a larger crystalline system. We assume that the occu- 
pied states can be represented in terms of well-localized 
orthonormal orbitals \cj>i), which we will refer to as Wan- 
nier functions. If wc introduce the velocity operator as 



- S [r,*], 



(1) 



then the total magnetic moment of the finite system in- 
volves the matrix elements (ipi\r x v^), where the 
are the occupied eigenstates of H . By invariance of the 
trace, this can be written in terms of WFs as 
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(2) 



where — e is the electron charge. The magnetization 
M can then be defined as the magnetic moment m 
per unit volume. For large but finite samples, all \<pi) 
that are sufficiently far from the surface become expo- 
nentially similar to bulk WFs, which we will denote 
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FIG. 1: (a) Four unit cells of the Haldane model. Filled 
(open) circles denote sites with Eo — — 2 (+2). Arrows in- 
dicate sign of the phase ip for second-neighbor hopping, (b) 
Net currents — ev^ associated with WFs, plotted at their cen- 
ters fj, for a 10x10 sample with (p = tt/4. Currents decrease 
rapidly into the bulk, so that only surface currents are visible. 

as \wi). For the electric polarization, the transition 
(<fii\r\(pi) — > (wi\r\wi) describes the polarization of peri- 
odic systems correctly [7]. Thus, it is tempting to assume 
that the magnetization should be expressible in a similar 
way in terms of the circulation (wi\r x v\v)i) of bulk WFs. 

We set out to verify this hypothesis in the context of 
numerical tight-binding calculations. For simplicity, we 
chose the Haldane model [11], which is comprised of a 
honeycomb lattice with two tight-binding sites per cell 
with site energies ±£b, real first-neighbor hoppings t\, 
and complex second-neighbor hoppings tie v as shown 
in Fig. 1(a). For our tests, we have chosen a lattice con- 
stant equal to unity, Eq = ±2, ti = 1 and ti — 1/3 and 
allowed ip to vary [16] . We treat the upper band as empty 
and the lower band as occupied. The corresponding WFs 
were obtained by acting with the band projector on a set 
of S- functions located on the sites with Eo — — 2 and 
applying a subsequent symmetric orthonormalization. 

The results of our numerical calculations for a finite 
30x30 sample are depicted as the symbols in Fig. 2. 
First, we calculated the total magnetic moment accord- 
ing to Eq. (2) and divided by the total sample area A to 
obtain the magnetization M indicated as circles in Fig. 2. 
Next, we evaluated the contribution to Eq. (2) from a sin- 
gle WF deep in the bulk of the finite sample and divided 
by the unit-cell area Aq to obtain a "local circulation" 
magnetization Mlc plotted as triangles in Fig. 2. We 
expected these two quantities to agree with each other 
within some numerical tolerance. On the contrary, the 
results indicate no agreement whatsoever. 

This surprising result forced us to reconsider our en- 
tire line of argument, revealing a profound oversight. It is 
easily shown that each bulk band of an insulating crystal 
must carry no net current, even in the absence of time- 
reversal symmetry. This means that each bulk-like WF, 
such as the one that we chose from the deep interior of the 
sample, must carry no net current, as is easily confirmed 
in the numerical calculation. We had assumed that the 
WFs at the boundary of the finite sample would likewise 
carry no net current, but this assumption is incorrect. In 
fact, the WFs near the boundary do carry a net current, 
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FIG. 2: Numerical results for the magnetization (in Bohr 
magnetons per unit area) for the Haldane model. Symbols 
denote results from a finite 30x30 sample, while curves rep- 
resent results of a k-space calculation on a periodic system. 



and the total circulation associated with these net cur- 
rents provides just the needed contribution to resolve the 
discrepancy. 

This is illustrated in Fig. 1(b), where we have plotted 
the net current — evj = — e(<t>i\v\(f>i) located at the Wan- 
nier center Fj = ((f>i\r\<pi) for each WF in the sample. 
While confirming that the net currents in the deep inte- 
rior are exponentially small, the results reveal that the 
WFs near the surface do carry a substantial current that 
contributes extensively to the total magnetic moment of 
the sample. Dividing by the sample area, we obtain an 
"itinerant circulation" contribution M\c to the magne- 
tization that is plotted as squares in Fig. 2. A glance 
at the figure suggests, and numerical tests confirm, that 
M = Mlc + Mio within numerical precision [17]. 

It is not surprising that Mlc , corresponding to the lo- 
cal circulation of a bulk WF, is a bulk property of the 
insulator [9]. It is far less clear whether M\c is also a 
bulk property. To show that it is, we go back to Eq. (2). 
For simplicity, we restrict ourselves henceforth to the 
case of a two dimensional system with a single occupied 
band (the generalization to three dimensions is straight- 
forward, while the multiband treatment is more subtle). 
Then Eq. (2) can be rewritten as 

M = ^ [ (Mr - ij) x v|&) + r, x (&|v|&) j , 

* LC IC 

(3) 

where again 'LC and 'IC correspond to local and itiner- 
ant circulation contributions, respectively. We divide the 
finite sample into an "interior" and a "surface" region in 
such a way that the latter occupies a non-extensive frac- 
tion of the total sample area in the thermodynamic limit. 
We label WFs from the surface and interior regions as 
\4> s ) and |R), respectively, where R is a lattice vector. 

Next, we note that the WFs in the surface region make 
a negligible contribution to the LC term of Eq. (3) since 
they occupy a non-extensive fraction of the area in the 
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FIG. 3: Horizontal slice from a sample that extends indefi- 
nitely in the vertical direction but is otherwise similar to the 
one in Fig. 1. Vertical dashed lines delimit bulk and surface 
regions in which WFs are labeled by s and s' , respectively. 



thermodynamic limit. The LC term then becomes 
M lc = -^-£<R|(r-r R )xv|R) 



R 



2A c 



(0|r x v|0) . 



(4) 



Here we have used translational symmetry and the fact 
that (R|v|R)=0 (since bulk bands carry no net current). 
Eq. (4) shows that Mlc can be expressed simply in terms 
of the bulk WF |0) in the home unit cell. Turning to 
the IC term, the interior WFs now make no contribution 
(again because (R|v|R)=0) so that 



(5) 



where the sum runs over surface WFs only. 

We now concentrate on this IC contribution and con- 
sider a vertical strip of which one horizontal section 
is sketched in Fig. 3, and choose vertical boundaries 
(dashed lines) to discriminate between interior and sur- 
face regions. Focusing on the right edge, we use labels s 
and s' to label WFs in the interior (x < x ) and surface 
(x > x ) regions, respectively. The vertical macroscopic 
current flowing in the right surface is 



e 

Al 



s ,y > 



(6) 



where the primed sum is further limited to WFs whose 
centers are inside a vertical segment of length Al, equal 
to the vertical repeat unit. The current carried by the 
i-th WF can be written in terms of contributions from 
its neighbors as 



--<&|[r,fl]|0i> 



(7) 



where vua = (2/7i) Ioiry Hji has the interpretation of 
a current "donated from WF j to WF i." Since v/ji) — 
—V{ij) , the total current carried by any subset of WFs 
can be computed as the sum of all vyy for which i is 
inside and j is outside the subset. Applying this to the 
piece of surface region considered above, Eq. (6) becomes 



(s,s'),y 



(8) 



Setting the boundary deep enough below the surface 
to be in a bulk-like region and invoking the exponen- 
tial localization of the WFs and of derived quantities 
like v/jj\, we can identify \<p s ) and \4> s >) with the corre- 
sponding bulk WFs. Exploiting translational symmetry, 
V (R,R'> = V( ,R'-R>, Eq. (8) becomes 
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(9) 



where the sum is still restricted to a segment of height AL 
The number of terms in Eq. (9) having a given value of 
R' R is just (R' x - R x )Al/A if (R' x -R x )>0 and zero 
otherwise. With a change of summation index, Eq. (9) 
becomes 



/„ = — 



2A 



(10) 



R 



where the factor of 2 enters because the sum has been 
extended to all R. For a boundary fragment of arbitrary 
orientation, Eq. (10) generalizes to I a = ^pGapnp, 
where n is the unit normal to the boundary, and 

G^ = --^-^Im(R|r a |0)(0|^|R)^. (11) 

R 

The contribution of this itinerant current to the magnetic 
moment (l/2c) §(r x I) dl is easily seen to be related to 
the antisymmetric part [18] of G, so that 
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(12) 



Equations (11) and (12) constitute our first major result, 
showing that the itinerant circulation contribution to the 
orbital magnetization can indeed be expressed as a bulk 
property in terms of the bulk WFs alone. 

In the remainder of this Letter, we show that the two 
contributions Mlc and M\c can both be converted into 
k-space expressions that can be evaluated directly in the 
Bloch representation. The WFs are defined via 



IR) 



(2nf 



d 2 k, 



ik-(r-R) 



(13) 



where |uk) = e~ tk ' r V'k) is the cell-periodic part of the 
Bloch function |V>k}- Inserting Eq. (1) into Eq. (4) and 
using r x r = 0, it follows that 



M LC 



lm(0|r x Fr|0) . 



2A hc 

Defining H k = e~ lk r He lk r and using that 

r |R) = l ^ j d 2 k e *-(r- R ) \d k U k ) , 



(14) 



(15) 



Eq. (14) becomes 



M LC = -r- Im 

2hc 



d 2 k 



(<9kUk| x H k |9ku k ) , (16) 
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in agreement with Ref. [9]. 

In order to convert M\c to k-space, we note that the 
matrix elements appearing in Eq. (11) are given by 



(0|r|R) 
<0|if|R) 



An 
(2ir) 2 

(2^) 2 



d 2 k A k e- 4k R , (17) 
d 2 k E k e- lk R , (18) 



where the Berry connection A k = i (u k |<9 k |u k ) is real and 
£k is the band energy. After some algebra including an 
integration by parts, we find 



G a fi 



d 2 k 
(2ir) 2 



Eu dkp A^ 



Inserting in Eq. (12) gives 



Mic 



e r d 2 i 
~2hcJ (2tt 



(19) 



(20) 



where f2 k = 9 k x A k is the Berry curvature. 

Interestingly, both magnetization contributions (16) 
and (20) are individually gauge- invariant, i.e., insensi- 
tive to the choice of phases of the Bloch functions used 
to construct the WFs. This was shown for Eq. (16) in 
Ref. [9], and it follows immediately for Eq. (20) because 
f2(k) is a gauge-invariant quantity. Each contribution is 
also invariant with respect to a shift of the zero of the 
Hamiltonian; the effect of such a shift is proportional to 
J d 2 kflk = 2ttC, where C is the Chern number which 
has been assumed to vanish. 

Adding Mlc and M\c, placing both in a common form, 
and returning to three dimensions, we find that the total 
magnetization of the crystalline solid can be expressed as 



M 



2hc 



Im 



d 3 k 
(2tt) 3 



<cW| x + Ek) |cW) . (21) 



This is our principal result. We note that Eq. (21) is 
consistent with Eq. (11) of Ref. [10], thereby providing a 
direct, fully quantum derivation of a result inferred there 
on the basis of semiclassical arguments alone. 

To confirm the correctness of our formulation, we used 
discretized [9, 19] versions of Eqs. (16) and (20) to cal- 
culate Mlc Me, and M for the Haldane model using a 
300x300 k-point mesh. The results, drawn as the lines 
in Fig. 2, are entirely consistent with the results for finite 
samples. We also carried out numerical tests which con- 
firmed that the extrapolation of the results from N x N 
finite samples for N — 10, 20 and 30 to N — > oo shows 
very precise agreement with the k-space calculation on 
a dense mesh. We can thus be confident that the for- 
mal derivations are correct and that there is no longer 
any possibility that terms in the magnetization are being 
overlooked. We have carried out similar tests for other 
tight-binding models with reduced symmetry [17], with 
similar results. 



In conclusion, we have derived a formula for the or- 
bital magnetization of a crystalline system by working in 
the Wannier representation, and we have demonstrated 
its correctness via numerical tests. While limited to the 
case of a non-interacting zero- Chern- number insulator in 
a vanishing (or commensurate) magnetic field, our re- 
sult nonetheless represents significant progress towards a 
more general theory of orbital magnetization. The re- 
sulting formula is easily evaluated in the context of a 
k-space electronic-structure code. The generalization to 
the multiband case will be discussed in a forthcoming 
publication. It remains tantalizingly uncertain whether 
such a Wannier-based approach can also be generalized 
to handle insulators with non-zero Chern numbers, met- 
als, or arbitrary magnetic fields. 
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